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Earthquake occurrence in nature is thought to result from correlated elastic stresses, leading to 
clustering in space and time. We show that occurrence of major earthquakes in California correlates 
with time intervals when fluctuations in small earthquakes are suppressed relative to the long term 
average. We estimate a probability of less than 1% that this coincidence is due to random clustering. 
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INTRODUCTION 

It is widely accepted QEHSHHIlSHIliQ 

that the observed earthquake scaling laws indicate the 
existence of phenomena closely associated with prox- 
imity of the system to a critical point. More specifi- 
cally, it has been proposed that earthquake dynamics 
are associated either with a second order critical point 
or a mean field spinodal [H,E3 that can 
be understood as a line of critical points. Mean field the- 
ories of the Ginzburg-Landau type have been proposed 
0,0 to explain the phenomenology associated with 

scaling and nucleation processes of earthquakes, which 
would in turn imply that a Ginzburg criterion is appli- 
cable 12]. If mean field Ginzburg-Landau equations do 
describe earthquakes, the dynamics must be operating 
outside the critical region, and fluctuations are corre- 
spondingly reduced. 



To summarize our results 

We compare the performance of two probability mea- 
sures that define the locations of future earthquake occur- 
rence: the spatially coarse-grained seismic intensity and 
the intensity change. We show that an order parameter 
^fj (t) can be defined based on the performance of these 
probability measures on a Receiver Operating Character- 
istic (ROC) diagram and that a generalized Ginzburg 
criterion Q(t) can be established measuring the relative 
importance of fluctuations in ^i(t). We find that since 
1960, major earthquakes in California with magnitudes 
m > 6 tend to preferentially occur during intervals of 
time when Q(t) < 1, consistent with mean field dynam- 
ics. Currently in northern California, Q(t) < 1. 



tudes — 124°E and — 115°E, coarse-grained in time inter- 
vals of one day. Only events above a magnitude threshold 
m,T > 3 are used to ensure catalog completeness. Fig- 
ure^shows the event locations. We tile the region with a 
spatially coarse-grained mesh of N boxes, or pixels, hav- 
ing side length 0.1°, about 11 km at these latitudes, ap- 
proximately the rupture length of an m ~ 6 earthquake. 
The average intensity of activity J(x, to, t 2 ) is constructed 
by computing the number of earthquakes n(x, to, £2) m 
each coarse-grained box centered at x since records be- 
gan at time to = 1932 until a later time t2 that will be 
allowed to vary: 7(x,to,t2) = n(x, tojta)- We then re- 
gard P M = P M (x, t 0l t 2 ) = I(x,to,fo)/ J J(x, t ,f2)dx as a 
probability for the location of future events m > tot for 
times t > f 2 • Previous work 0, Q, indicates that 
P M is a good predictor of locations for future large events 
having m > 5. 

The intensity change map builds upon the intensity 
map by computing the average squared change in inten- 
sity over a time interval At = t 2 — t\. Here we use 
At = 13 years [T3L ITij ]. We compute n(x,t(,,ti) and 
n(x, tfejta) for the two times t\ and t2, where t2 > ti, 
beginning at a base time t&, where ti > % > to. Comput- 
ing the change in numbers of events as An(x, tb, ti , t2) = 
n(x, tfe,t2) — n(x,tb,ti), we then define the intensity 
change AJ(x, t\,t2) by normalizing An(x, tb, ti, t2) to 
have spatial mean zero and unit variance, yielding 
An'(x, tb, ti, t2), and then averaging An'(x, t&, ti, t2) 
over all values for tb from to to t\: A7(x, ti,t2) =< 
An'(x, tb, fi, f2) >t b - The corresponding probability is 
P A = P A (x,t 1 ,t^^AJ(x,t 1 ,t2)] 2 / J[A/(x,t 1 ,t 2 )] 2 dx. 
Previous work [13|, [lj, [lj| has found that Pa is also a 
good predictor of locations for future large events having 
to > 5. Pa can be viewed as a probability based upon 
the squared change in intensity. 



INTENSITY MAPS AND INTENSITY CHANGE 
MAPS 

The data set we use is the ANSS catalog of earthquakes 
[Lsj ] between latitude 32°N and 40°N and between longi- 



BINARY FORECASTS 

Binary forecasts are a well-known method for con- 
structing forecasts of future event locations and have 
been widely used in tornado and severe storm forecast- 
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FIG. 1: Map of earthquake (m > 3) epicenters in California 
from 1932 to the present. Circles are events with m > 6 
since 1960. Red epicenters define the area used to analyze 
seismicity in northern California; blue epicenters define the 
area used for southern California. 
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FIG. 2: (A) Pa (intensity change) map (left) and P M (average 
intensity) map (right). A decision threshold D was chosen 
leading to 75 hotspots in each map. Maps were computed for 
t =1 December 1994 with t^—ti = 13 years. (B) ROC curves 
for Pa and P M corresponding to Figure Q Here we have used 
m c = 4. 



ing [lot llq . We construct binary forecasts for m > m c 
and for times t > t 2 , where m c is a cutoff magnitude. In 
past work (l3lll4lllq we have taken m c = 5, but we now 
remove this restriction. In our application, the prob- 
abilities P M = P M (x, to,t 2 ) and Pa = Pa(x, are 
converted to binary forecasts P M = B^(D, x, to, t 2 ) and 
Pa = Ba{D,x, tijta) by the use of a decision threshold 
-D, where D £ [0,max{P M }] or D 6 [0,max{PA}] respec- 
tively mm 



For a given value of D, we set P M = 1 where P^ > D 
and P,, = otherwise. Similarly, we set Pa = 1 where 
Pa > D and Pa = otherwise. The set of pixels 
{x M (D)} where P M = 1 and {xa(-D)} where Pa = 1 
then constitute locations where future events m > m c are 
considered to be likely to occur. We call these locations 
hotspots. The locations where B^ = and Pa = are 
sites where future events m > m c are unlikely to occur. 
In previous work, intensity maps and intensity change 
maps at a particular value of D were called Relative In- 
tensity maps and Pattern Informatics maps. Examples 
of binary forecast maps are shown in Figure 



RECEIVER OPERATING CHARACTERISTIC 
(ROC) DIAGRAMS 

A series of P>-dependent contingency tables are con- 
structed using the set of locations {x g (m c )} where the 
q = 1 , . . . , Q large events m > m c are observed to actu- 
ally occur during the forecast verification period t > t 2 . 
The contingency table has 4 entries, a — ► d, whose values 
are determined by some specified rule set |l5l Il6| . Here 
we use the following rules for given D (same rules for 
both and "A" subscripts): 

1. a is the number of boxes in {x(P/)} which are also 
in {x g (m c )} 

2. b is the number of boxes in {x(D)} whose location 
is not in {x (J (m c )}, i.e., is in the complement to 
{x q (m c )} 

3. c is the number of boxes in the complement to 
{x(P>)} whose location is in {x g (m c )} 

4. d is the number of boxes in the complement to 
{x(P/)} whose locations are in the complement to 
{x,(m c )} 
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The hit rate is then defined as H = a /(a + c), and the 
false alarm rate is defined as F = 6/(6 + d). Note that 
with these definitions, a + c = Q, a + b = number of 
hotspots, and a + b + c + d = N. 

The ROC diagram 0, 0] is a plot of the points 
{H, F} as D is varied. Examples of ROC curves cor- 
responding to the intensity and intensity change maps 
in Figure are shown in Figure 03. A perfect forecast 
of occurrence (perfect order, no fluctuations) would con- 
sist of two line segments, the first connecting the points 
(H, F) = (0,0) to (H,F) = (1,0), and the second con- 
necting (H,F) = (1,0) to (H,F) = (1,1). A curve of this 
type can be described as maximum possible hits (H = 1) 
with minimum possible false alarms (F — 0). Another 
type of perfect forecast (perfect order, no fluctuations) 
consists of two lines connecting the points (0, 0) to (0, 1) 
and (0, 1) to (1, 1), a perfect forecast of non- occurrence. 

The line H — F occupies a special status, and corre- 
sponds to a completely random forecast [1.4 Il6| (max- 
imum disorder, maximum fluctuations) where the false 
alarm rate is the same as the hit rate and no information 
is produced by the forecast. Alternatively, we can say 
that the marginal utility |l7| of an additional hotspot, 
dH/dF, equals unity for a random forecast. 

For a given time-dependent forecast H(F,t), we con- 
sider the time-dependent Pierce Skill Score H(F, t) — F 
, which measures the improvement in performance of 
H(F, t) relative to the random forecast H = F. A Pierce 
function ^(t) measures the area between H(F, t) and the 
random forecast: 



*(*) 



(H(F, t) - F)dF 



(1) 



= A(t) - F^J2, 



where 



A(t) 



H(F,t)dF. 



(2) 



The upper limit F max on the range if integration is a 
parameter whose value is set by the requirement that the 
marginal utility 01 of the forecast of occurrence H(F, t) 
exceeds that of the random forecast H — F: 



d 

dF 



{H(F,t) - F} > 0. 



(3) 



Since H(F, t) curves are monotonically increasing, F max 
is determined as the value of F for which dH(F, t)/dF = 
1. For the forecasts we consider, we find that F max w 0.2, 
as can be seen from the examples in Figure |2J3. 



ORDER PARAMETER AND GENERALIZED 
GINZBURG CRITERION 

We define an order parameter as the Pierce func- 
tion ^ T (t) obtained using as the probability P T = 



P T (x, <i,< 2 ) = n(x,ti,t 2 )/ J n(x, ti,i 2 )dx, where P T is 
the average normalized intensity of seismic activity dur- 
ing t\ to t 2 . Using P T and the decision threshold D, we 
construct a binary forecast B T = B T (D, x, t\, i 2 ) • Eval- 
uating the forecast B T during the time interval i 2 to t 
produces the ROC diagram H T (F,t). For the case of 
forecasts having positive marginal utility relative to the 
random forecast, ^ T (t) > 0. If past seismic activity is un- 
corrected with future seismic activity, P T is equivalent 
to a random forecast, and *£ T (t) = 

Corresponding to the order parameter (t), we define 
a function Q(t) to indicate the relative importance of fluc- 
tuations with respect to forecasts of occurrence. We note 
that the probability V^a is a measure of the mean squared 
change of intensity, a measure of fluctuations in seismic 
intensity, during t\ to i 2 , and that the probability P^ is 
a measure of the average intensity over the entire time 
history (to to t 2 ). We will refer to Pa as the "fluctuation 
map" or "change map" , and P M as the "average map" . 

Using the corresponding ROC functions we define 



(4) 



where \& a(£) is based upon the ROC curve computed 
using Pa, {Ha{F, t),F} and ^(t) is based upon the 
ROC curve computed using P^, {H^(F, t), F}. We can 
say that when Q(t) < 1, "fluctuations are less significant 
relative to the mean" in the sense that the fluctuation 
map provides a poorer forecast than the mean map. This 
statement is equivalent to the Pierce difference function: 



AA(t) = A^t) - A A (t) >0. 



(5) 



This difference function can be considered to be a gener- 
alized Ginzburg criterion [l^. 

To examine these ideas, we compare a plot of Q(t) 
with activity of major earthquakes (m > 6) in Califor- 
nia. We first consider the Gutenberg-Richter frequency- 
magnitude relation / = 10° • 10~ hm , where / is the num- 
ber of events per unit time with magnitude larger than m 
and a and b are constants, a specifies the level of activity 
in the region, and 6=1. 

To construct ROC curves, we consider t to be the cur- 
rent time at each time step and test the average map 
and change map by forecasting locations of earthquakes 
during £ 2 to t. We use events having m > mr, where 
rriT is some threshold magnitude. Note that / _1 speci- 
fies a time scale for events larger than m: 1 event with 
m > 6.0 is associated on average with 10 m > 5.0 events, 
100 m > 4.0 events, etc. Without prior knowledge of 
the optimal value for tut, we average the results for a 
scale-invariant distribution of 1000 tot > 3.0 events, 794 
tot > 3.1 events, 631 tot > 3.2 events, . . ., 10 tot > 5.0 
events. We terminate the sequence at tot > 5.0 due 
to increasingly poor statistics. To control the number 
of earthquakes with to > tot in the snapshot window 
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(t2 to t), we determine the value of ti that most closely 
produces the desired number of events within the snap- 
shot window. It is possible to have fluctuations in actual 
number of events if the snapshot window includes the oc- 
currence time of a major earthquake, when there may be 
many events to > tot in the coarse-grained time intervals 
of length 1 day following the earthquake. 

A central idea is that the length of the snapshot win- 
dow is not fixed in time; it is instead fixed by earthquake 
number at each threshold magnitude = 3.0,3.1,3.2, 
and so forth. Nature appears to measure "earthquake 
time" in numbers of events, rather than in years. "Earth- 
quake time" is evidently based on stress accumulation 
and release, that is, earthquake numbers, rather than in 
months or years jlOj . 

Results are shown in Figure for the region of Cal- 
ifornia shown in Figure Q At top of either plot is the 
Pierce difference function AA(t) = A^it) —A^(t), and at 
bottom is earthquake magnitude plotted as a function of 
time from 1 January 1960 to 31 March 2006. The vertical 
lines in each top panel are the times of all events to > 6 in 
the region during that time interval. It can be seen from 
Figures n] an d El that there are 11 to > 6 events in north- 
ern California and 10 such events in southern California. 
For both areas, these major events are concentrated into 
8 distinct "episodes" corresponding to 8 main shocks. In 
each plot, 7 of the 8 major episodes fall during ("black") 
time intervals where AA(t) > 0, or they either begin or 
terminate such a time interval. If a binomial probability 
distribution is assumed, the chance that random cluster- 
ing of these major earthquake episodes could produce this 
temporal concordance can be computed. For Figure EK, 
where black time intervals constitute 36.8% of the to- 
tal, we compute a 0.46% chance that the concordance is 
due to random clustering. For Figure EP, the respective 
numbers are 19% of the total time interval, and 0.0058% 
chance due to random clustering. Our results support the 
prediction that major earthquake episodes preferentially 
occur during time intervals when fluctuations in seismic 
intensity, as measured by ROC curves, are less important 
than the average seismic intensity. 
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FIG. 3: Value of the Pierce difference function AA(t) (top) 
and magnitude (bottom) as a function of time for events oc- 
curring on the map area of Figure^ Vertical black lines repre- 
sent times of major earthquakes having m > 6 in the respec- 
tive regions. Differences are computed for a scale-invariant 
distribution of magnitude thresholds in the snapshot window 
from mi = 3.0 to mi = 5.0. Area integration is performed 
for F £ [0.0,0.2]. (A) Northern California (red epicenters in 
Figure 0. (B) Southern California (blue epicenters in Fig- 
ure 0. 
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